function [e_est, h_est, Delta_Bound_mod] = bound_demandcomp(y_start, p_start, p_final, points, g, g2, g22)

delta = (p_final - p_start)/(points-1);

e_est = NaN*ones(points,1); 
e_est(1,1) = y_start;
h_est = NaN*ones(points,1); 
h_est(1,1) = g(1,1);

for ii = 1:points-1
    e_est(ii+1,1) = e_est(ii,1) + delta * (g(ii,1) + (e_est(ii,1)-y_start)*g2(ii,1));
    h_est(ii+1,1) = g(ii+1,1) + (e_est(ii+1,1)-y_start)*g2(ii+1,1);
end

Delta_Bound = NaN*ones(points,1); Delta_Bound(1,1) = 0;
Delta_Bound_mod = NaN*ones(points,1); Delta_Bound_mod(1,1) = 0;

for ii = 1:points-1
    Delta_Bound (ii+1,1) = (1 + delta * g2(ii,1))*Delta_Bound(ii,1) - 0.5*delta*(e_est(ii,1)-Delta_Bound(ii,1)-y_start)^2 * g22(ii,1);

    term1 = (e_est(ii,1)-Delta_Bound_mod(ii,1)-y_start)^2;
    term2 = (e_est(ii,1)-y_start)^2;
    Delta_Bound_mod (ii+1,1) = (1 + delta * g2(ii,1))*Delta_Bound_mod(ii,1) - 0.5*delta*(max(term1,term2)) * g22(ii,1);
end
%
